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ABSTRACT 


The  exergy  method  of  systems  analysis  uses  both  the  First  and  Second  laws  of 
Thermodynamics  to  determine  where  losses  occur.  This  method  has  been  shown  to  be  superior  to 
conventional  heat  balance  analyzes.  A  primary  contributor  to  the  exergy  decrease  is  the  production 
of  entropy.  Entropy  is  produced  in  shock  waves  and  boundary  layers,  as  well  as  in  other  regions  of 
the  flowfield.  A  method  of  analysis  is  presented  to  quantify  the  entropy  that  is  produced  in  the  bow 
shock  wave  and  the  boundary  layer  of  a  space  launch  vehicle  for  Mach  numbers  between  1 . 1  and  10. 
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I. 


INTRODUCTION 


As  commercial  interest  in  space  launch  vehicles  has 
increased  over  the  last  twenty  years,  there  has  been  an 
increased  emphasis  on  improving  the  efficiency  of  vehicles 
used  for  space  launch  services.  This  emphasis  will  continue 
to  receive  attention  as  competition  forces  the  providers  of 
launch  services  to  decrease  costs. 

This  thesis  is  an  initial  step  in  the  development  of  an 
exergy  methodology  tool  for  the  conceptual,  preliminary,  and 
detailed  design  of  aerospace  vehicles. 

In  order  to  determine  the  efficiency  of  a  particular 
launch  system,  it  is  necessary  to  identify  the  various 
sources  of  efficiency  loss  and  the  quantity  of  efficiency 
loss  associated  with  each  source  within  the  system.  Although 
energy  cannot  be  created  or  destroyed,  it  can  be  transformed 
into  unavailable  forms.  This  thesis  will  present  an 
investigation  of  one  method  that  may  be  used  to  examine  the 
efficiency  loss  from  two  specific  areas  of  a  typical  launch 
vehicle.  This  method  uses  the  concept  of  available  energy  or 
exergy  as  the  basis  for  determining  the  decrease  in  the 
energy  state  that  occurs  in  the  shock  wave  and  the  boundary 
layer  of  a  supersonic  vehicle. 

Typical  performance  analyzes  of  energy  systems  use  only 
the  first  law  of  thermodynamics.  This  method  of  analysis  is 
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based  on  the  principle  of  conservation  of  energy,  and 
suggests  that  an  energy  balance  accounting  for  all  energy 
flow  into  or  out  of  the  system  will  provide  a  measure  of  the 
system  efficiency.  While  the  first  law  approach  is  extremely 
important  in  the  performance  analysis  of  energy  systems,  it 
provides  only  the  overall  performance  or  the  actual  gross 
performance  of  the  system  or  process.  To  establish  the 
magnitude,  location,  and  type  of  loss  within  a  system  it  is 
necessary  to  use  the  Second  law  of  Thermodynamics.  The 
second  law  accounts  for  the  irreversibility  of  all  real 
processes  and  thus  provides  a  method  of  calculating  the 
losses  associated  with  entropy  production  that  occurs. 

The  exergy  method  of  analysis  uses  both  the  First  and 
Second  laws  of  Thermodynamics  to  determine  the  energy  that 
is  available  at  various  points  in  a  system.  While  the  exergy 
method  of  analysis  has  been  applied  to  conventional  closed 
cycle  thermal  systems  such  as  power  generation  and 
refrigeration  systems,  it  has  not,  to  the  author's 
knowledge,  been  applied  to  the  open  cycle  space  launch 
vehicle . 

"Exergy  is  defined  as  the  work  that  is  available  in  a 
gas,  fluid,  or  mass  as  a  result  of  its  nonequilibrium 
condition  relative  to  some  reference  condition. " [Ref .  l:p. 
33]  Since  exergy  is  defined  as  the  available  work  in  a  gas, 
fluid,  or  mass,  it  is  a  special  case  of  defining  available 
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energy.  Exergy  is  an  explicit  property,  and  therefore  can  be 
calculated  at  any  point  in  an  energy  system. 

The  exergy  method  permits  the  examination  of  any  process 
in  relation  to  the  theoretically  most  efficient  manner  that 
the  process  could  be  carried  out  within  the  environment.  An 
exergy  analysis  allows  one  to  quantify  the  loss  of 
efficiency  in  a  process  due  to  the  loss  of  the  quality  of 
the  energy. 

Exergy  is  calculated  relative  to  a  reference  condition 
by  the  following  general  equation: [Ref .  l:p.  34] 

Exergy  =  <u-u0)  -  T0(s-s0)  +  lS(v-va)  *  * 

C 

<z~zo)  yr~T  +  E  +  (!) 

Qc'*  o 

EiAiFi  (3T*-tI~4T0T2)  +  ... 

where  the  subscript  o  denotes  the  reference  condition.  The 
terms  on  the  right  side  of  the  equation  are  internal  energy, 
change  in  heat  energy,  flow  work,  kinetic  energy  per  unit 
mass,  potential  energy,  chemical  energy,  and  radiation 
energy  emission  respectively.  In  the  analysis  of  specific 
systems  or  processes  one  or  more  of  the  above  terms  may  be 
omitted,  and  only  those  relevant  to  the  problem  being 
considered  need  be  included. [Ref  l:pp. 34-36] 

This  thesis  focuses  on  the  quantity  of  exergy  decrease  as 
a  result  of  entropy  production  across  the  shock  wave  and  as 
a  result  of  boundary  layer  formation.  As  can  be  seen  from 
equation  1,  an  increase  in  entropy  results  in  a  decrease  in 
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exergy.  The  Second  law  of  Thermodynamics  tells  us  that  all 
real  process  are  irreversible  and  that  entropy  is  increased 
in  an  irreversible  process. 
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II.  EXERGY  DECREASE  ACROSS  THE  SHOCK  WAVE 


Supersonic  flow  over  a  blunt  body  will  result  in  the  bow 
shock  wave  being  detached  from  the  body  as  shown  in  Figure 
1.  The  shape  of  the  detached  shock  wave,  its  detachment 
distance,  and  the  flowfield  between  the  body  and  the  shock 
are  determined  by  the  shape  of  the  body  and  the  free  stream 
Mach  number . [Ref . 2 : p . 128 - 129 ] 
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Since  the  velocity  of  the  body  is  supersonic,  the  flow  ahead 
of  the  body  is  unaware  of  the  presence  of  the  body  until  the 
shock  wave  is  encountered.  The  shock  wave  directly  in  front 
of  the  body  is  approximately  normal  to  the  free  stream  and 
consequently  the  flow  behind  the  shock  wave  at  this  point  is 
subsonic.  The  subsonic  flow  region  is  defined  by  the  shock 
wave,  the  body,  and  the  sonic  line  as  shown  in  Figure  l. 

Because  the  flowfield  between  a  blunt  body  and  its 
detached  shock  wave  contains  both  subsonic  and  supersonic 
flow,  an  exact  determination  of  the  flowfield  dynamics  is 
extremely  complicated  to  formulate.  The  steady  state 
flowfield  in  the  subsonic  region  is  described  by  elliptic 
equations  ard  in  the  supersonic  region  by  hyperbolic 
equations.  The  sonic  line  divides  these  two  regions.  It  is 
the  problem  of  predicting  the  flowfield  in  the  transition 
region  between  the  subsonic  and  supersonic  flow  that  makes 
the  flowfield  calculations  for  a  supersonic  blunt  body  so 
difficult . 

The  flowfield  behind  a  detached  shock  wave  cannot  be 
precisely  determined  without  the  application  of 
sophisticated  numerical  methods.  In  keeping  with  the 
development  of  a  design  tool,  a  simplified  method  for 
determining  the  approximate  shape  of  the  detached  shock  wave 
that  was  developed  by  Moeckel  [Ref.  3]  is  used.  This  method 
was  subsequently  examined  by  Love  [Ref.  4]  and  with  some 
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modification  was  found  to  be  valid  for  blunt  bodies 
traveling  at  supersonic  and  hypersonic  speeds. 

Love  introduced  two  major  modifications  to  Moeckel's 
method.  The  first  of  these  was  in  the  method  used  to 
calculate  the  shock  detachment  distance  from  the  body. 
Moeckel  used  a  continuity  relation  to  calculate  detachment 
distance  and  shock  shape.  Love  calculates  shock  detachment 
distance  by  the  following  equation: 


x 


/ 


0.5  C  cot  $det 


(2) 


where  x'/d'  is  the  nondimens ional  detachment  distance,  C  is 
a  factor  that  depends  on  type  of  body  (spherical, 
cylindrical,  etc.),  and  5det  is  the  required  deflection 
angle  for  shock  wave  detachment,  from  a  conical  body,  for 
the  given  Mach  number  (see  Figure  .2)  . 

The  second  modification  made  by  Love  was  in  the 
calculation  of  the  angle  eta  (tj)  that  the  sonic  line  makes 
with  the  freestream.  Love  conducted  wind  tunnel  tests  and 
compared  his  experimental  results  with  his  theoretical 
calculations  [Ref .4:pp. 13-15] .  For  this  analysis,  the  angle 
r]  corresponding  to  l.lsMslO  was  obtained  by  developing  the 
followirg  analytical  expression  for  the  values  of  r\  vs  Mach 
number  as  set  forth  by  Love  [Ref.4:p.45] : 

,  3.5901  +  56.9298  *log(M) 

1  ~  57.29578 
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Figure  2  s  Representation  of  Flow  With  Detached  Shock 
and  Notation  Used  in  Analysis. 

A.  BODY  SHAPE 

The  shape  of  the  deflecting  body  is  one  of  the  factors 
determining  the  shock  wave  shape.  For  this  study,  the  body 
was  patterned  somewhat  loosely  after  the  Titan  III -A  launch 
vehicle.  This  body  shape  is  essentially  a  blunt  cone  and 
cylinder  arrangement.  The  dimensions  of  the  launch  vehicle 
used  in  the  analysis  are  shown  in  Table  1. 


TABLE  1:  LAUNCH  VEHICLE  DIMENSIONS 


Length  Overall 

120  Ft. 

Diameter 

10  Ft. 

Nose  Radius 

1.5  Ft. 

Nose  Cone  Length 

10  Ft. 
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B .  MACH  NUMBER 


Besides  the  shape  of  the  body,  the  other  factor  in 

determining  the  shape  of  the  shock  wave  is  the  freestream 

Mach  number.  This  study  determined  the  shock  shape  based  on 
l.lsMslO.  Although  this  is  far  short  of  the  Mach  number 
required  to  reach  orbit,  the  purpose  of  the  study  was  to 

demonstrate  the  validity  of  the  method,  which  can  be 

accomplished  in  the  listed  Mach  range. 

C.  ASSUMPTIONS 

The  method  of  approximating  the  shock  shape  and  location 
developed  by  Moeckel  is  based  on  two  assumptions: 

(1)  The  form  of  the  shock  between  its  foremost  point 
and  its  sonic  point  is  adequately  represented  by 
a  hyperbola  asymptotic  to  the  free  stream  Mach 
lines;  and 

(2)  the  sonic  line  between  the  shock  and  body  is 
straight  and  inclined  at  an  angle  that  depends 
only  on  the  freestream  Mach  number.  [Ref.3:p.l] 

The  first  assumption  states  that  the  shock  shape  between  its 

vertex,  or  foremost  point,  and  its  sonic  point  may  be 

represented  by  a  hyperbola.  For  this  study,  the  assumption 

was  made  that  this  approximation  for  the  shock  shape  is 

valid  past  the  sonic  point  of  the  shock  wave.  For 

calculation  purposes,  the  shock  wave  was  taken  from  the 

vertex  of  the  wave  to  the  point  where  the  angle  between  the 

shock  wave  and  the  asymptote  that  defines  the  hyperbola 

differ  by  no  more  than  one  per  cent.  The  equation  for  the 

hyperbola  used  in  these  assumptions  is 
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(4) 


where  x  is  the  longitudinal  axis  of  the  body,  y  is  the  axis 
perpendicular  to  the  x  axis,  0  is  the  cotangent  of  the  Mach 
angle,  and  xQ  is  the  distance  from  the  vertex  of  the  shock 
wave  to  the  intersection  of  its  asymptotes . [Ref . 3 :p. 1] 

The  second  assumption  made  by  Moeckel  concerning  the  sonic 
line  was  modified  by  Love  as  previously  discussed  and  was 
used  in  this  analyzes. 

Air  is  the  fluid  used  for  all  of  the  analysis  presented 
in  this  thesis.  However,  the  analysis  method  is  equally 
valid  for  any  gas,  and  would  only  require  using  the 
appropriate  value  for  the  ratio  of  specific  heats  (y)  and 
the  new  specific  gas  constant.  It  has  been  assumed  that  the 
fluid  behaves  as  a  perfect  gas  in  all  of  the  calculations. 
This  assumption  is  valid  for  velocities  up  to  approximately 
M=3  to  M=5  [Ref  5:pp. 17-20].  Above  these  velocities 
dissociation  of  the  fluid  begins  to  occur  due  to  atmospheric 
heating.  If  the  gas  continues  to  be  heated,  ionization  will 
occur.  As  the  degree  of  dissociation  and  ionization 
increases,  the  perfect  gas  assumption  becomes  less  valid. 

The  calculation  of  total  entropy  across  the  shock  wave 
requires  a  determination  of  the  mass  flow  rate  of  air 
through  the  shock  wave.  The  air  mass  flow  rate  calculation 
is  based  on  the  density,  cross  sectional  area  being  studied, 
and  the  velocity  of  the  air.  Since  the  density  of  air  in  the 
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atmosphere  varies  with  height  above  sea  level,  an  arbitrary 
height  of  30,000  meters  was  assumed  for  all  M  values.  The 
density  and  temperature  of  air  at  this  altitude  was  taken 
from  the  International  Civil  Aviation  Organization  (ICAO) 
standard  atmosphere  tables  [Ref  6] .  The  calculations  may  be 
used  for  any  altitude  by  inserting  the  appropriate  density 
and  temperature . 

D .  EXERGY  CALCULATION 

From  equation  (1) ,  and  neglecting  all  terms  on  the  right 
except  the  entropy  term,  the  total  exergy  decrease  through 
the  shock  wave  is  proportional  to  the  increase  in  entropy 
across  the  shock.  The  entropy  change  across  a  shock  is  given 
by  [Ref.  7 : p .  10} 

As  2TMfslnJ«-W-l)  , 

cv  7  +  l  (y  - 1)  sin2  6+2  {  1 

where  cv  is  the  specific  heat  at  constant  volume  for  air,  y 
is  the  ratio  of  specific  heats,  M  is  the  freestream  Mach 
number,  and  6  is  the  local  shock  wave  angle  with  reference 
to  the  freestream  direction.  By  using  the  method  of  Moeckel, 
modified  by  Love,  to  determine  the  approximate  shock  shape, 
the  entropy  increase  across  the  total  shock  was  calculated 
by  equation  5  and  numerical  integration  using  cumulative 
applications  of  Simpson's  Rule.  The  number  of  integrations 
varied  with  the  shock  strength  and  shape  running  from  a  low 
of  50  for  Mach* 1.1  to  a  high  of  2513  for  Mach=10. 
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The  entropy  calculations  were  performed  by  the  Fortran 
77  program  contained  in  Appendix  A.  For  the  assumed  body 
shape,  this  program  calculates  the  approximate  shock  wave 
shape  and  the  total  rate  of  entropy  change  across  the  shock 
for  a  designated  control  volume  and  designated  Mach  numbers 
between  1.1  and  10.  The  calculations  are  carried  to  three 
decimal  places  with  the  knowledge  that  assumptions  and 
round-off  error  may  have  affected  the  accuracy.  The  results 
are  summarized  in  Table  2. 


TABLE  2  s  ENTROPY  CHANGE  ACROSS  SHOCK 


Mach  Number 

Rate  of  Entropy  Change 
Across  Shock  (KJ/sec} 

l.l 

.  .  9 

2.0 

50.803 

3.0 

301.469 

4.0 

836.941 

5.0 

1662.150 

6.0 

2764.771 

7.0 

4131.494 

8.0 

5759.238 

9.0 

7641.229 

10.0 

9774.107 

As  can  be  seen  from  equation  5,  the  entropy  rise  across 
the  shock  increases  as  the  shock  wave  angle  (0)  ,  and  the 
freest ream  Mach  number  ( M x)  increase.  For  a  detached  shock 
wave,  the  shock  wave  angle  is  greatest  immediately  ahead  of 
the  body  where  it  is  approximately  normal  to  the  freestream 
and  decreases  as  it  curves  around  the  body  until  it 
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eventually  asymptotes  the  freestream.  This  indicates  that 
for  a  given  body  shape  and  Mach  number,  the  greatest  exergy 
loss  occurs  in  the  vicinity  of  the  vertex  of  the  shock  wave. 
This  is  shown  graphically  in  Figure  3.  The  vertical  axis  is 
shown  as  entropy  change  over  the  specific  gas  constant  (R' ) , 
and  the  horizontal  axis  is  shown  as  the  angle  phi  ( <f> )  ,  which 
is  90  -  6. 


Figure  3:  Entropy  Rise  Across  Shock  Wave. 


E.  VERIFICATION  OF  THE  COMPUTER  MODEL 

Comparisons  between  theoretical  shock  location  using  the 
program  of  Appendix  A  and  actual  shock  location  as 
determined  by  Love  for  Mach  numbers  of  1.94  and  6.8  for  a 
combination  sphere  and  cylinder  body  shape  are  shown  in 
Figures  4  and  5  respectively.  These  comparisons  show  good 


13 


agreement  between  calculated  shock  shape  and  shock  location. 
The  data  points  corresponding  to  experimental  results  were 
taken  from  charts  [Ref.  4:pp  46,51]  and  therefore,  precise 
quantitative  evaluations  regarding  goodness  of  fit  are  not 
possible.  However,  the  experimental  data  appear  to  be  within 
five  percent  of  the  predicted  values  of  the  shock  profile. 


Shock  Shape.  Mach  1.94 


Shock  Shape.  Mach  6.8 
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Another  reality  check  for  the  computer  model  was  made  by 
calculating  the  mass  flow  of  air  through  the  control  volume 
defined  by  the  shock  wave  at  l.lsMslO  at  an  altitude  of 
30,000  meters.  Since  the  vehicle  is  supersonic,  there  is  no 
deflection  of  the  flow  around  the  body  prior  to  passing 
through  the  shock.  The  validity  of  the  shock  shape  generated 
by  the  program  is  good  as  shown  in  Figures  4  and  5  above. 

The  shock  wave  may  be  thought  of  as  a  blunt  cone.  If  the 
cone  is  terminated  at  any  point,  the  distance  from  the 
centerline  to  the  shock  is  the  radius  of  the  base.  The  mass 
flow  through  the  area  of  the  base  can  be  easily  calculated. 
This  was  done  and  the  result  was  compared  with  the  mass  flow 
calculated  by  the  program  using  the  same  Simpson's  rule 
integration  technique  that  was  used  to  calculate  the  entropy 
rise  across  the  shock.  The  results  of  these  two  calculations 
differed  by  less  than  three  per  cent  for  all  Mach  numbers. 
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III.  EXERGY  DECREASE  IN  THE  BOUNDARY  LAYER 


If  a  control  volume  is  taken  around  a  cylindrical  body 
moving  through  a  gas,  there  will  be  a  change  in  the  momentum 
of  the  gas  taken  upstream  of  the  body  from  that  taken 
downstream  of  the  body.  This  difference  in  momentum  is  a 
function  of  the  drag  on  the  body,  and  is  given  by  [Ref. 

8 : pp .  23-24] 


D  =  /  pu(U-u)  dA  (6) 

J  A 

where  the  integral  is  taken  over  the  end  surface  of  the 
cylinder  at  a  large  distance  downstream,  U  is  the  velocity 
upstream,  and  u  is  the  velocity  at  an  arbitrary  point  in  the 
wake.  Letting 


U-u 


According  to  Crocco's  theorem  ^ef.  9] 


(7) 


ay  dz 


T 

u 


( —  c?y- 
dy 


ds 

dz 


dz) 


(8) 


Substituting  equations  7  and  8  into  equation  6,  we  get 


D=  f  p udAf  —  (-f^  dy+  dz)  (9) 

j  a  Jo  u  ay  dz 

Replacing  pudA  {the  mass  flow  through  the  cross -sectional 
area  dA)  by  dm  and  T/u  by  Tm/U,  which  is  justified  if  second 
order  terms  are  neglected  [Ref.  8:p.24],  in  equation  9  gives 


16 


D  =  —fdmls-sj 
U  J  A 


(10) 


Taking  the  integral  over  the  area  of  the  wake,  equation  10 
can  be  written  as 

DU  =  Tm — —  (11) 

sec 

This  is  an  application  of  the  theory  of  Oswatitsch  [Ref. 

10]  .  The  freestream  temperature  and  velocity  are  known  or 
can  be  calculated  using  the  ICAO  tables  [Ref.  6]  for  a  given 
altitude  and  Mach  number.  If  the  drag  can  be  calculated, 
then  the  entropy  increase  can  be  calculated  using  equation 
11.  To  calculate  the  profile  drag  of  the  body,  the  flow 
parameters  at  the  boundary  layer  edge  and  the  boundary  layer 
shape  are  needed. 

The  theory  that  the  flow  surrounding  an  immersed  body 
could  be  divided  into  two  parts  was  first  proposed  by  L. 
Prandtl  in  1904.  Prandtl  stated  that  the  viscosity  in  the 
major  part  of  the  flow  field  has  a  negligible  effect  on  its 
motion  and  consequently  in  most  cases  could  be  ignored.  The 
other  part  of  the  flow  field  consisted  of  a  narrow  region 
next  to  the  body  where  the  viscous  forces  are  not 
negligible.  [Ref. 11] 

This  narrow  region  near  the  body  is  known  as  the 
boundary  layer.  In  the  boundary  layer  the  velocity  of  the 
fluid  increases  rapidly  from  zero  at  the  surface  to  a  value 
which  then  changes  slowly  with  increasing  distance  from  the 
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body.  AC  an  infinite  distance  from  the  body,  the  boundary 
layer  and  the  freestream  merge  smoothly  together.  For 
practical  applications,  the  boundary  layer  conditions 
approach  the  freestream  condition  very  rapidly  within  a 
relatively  small  distance  from  the  surface  as  shown  in 
Figure  5. 


If  the  freestream  velocity  is  taken  as  U,  then  the  edge  of 
the  boundary  layer  may  be  designated  as  that  distance  from 
the  body  where  the  velocity  of  the  stream  is  .995U  or  ue. 

Boundary  layers  may  be  either  laminar  or  turbulent.  In 
laminar  flow,  the  flow  within  the  boundary  layer  is 
characterized  by  smooth  streamlines  approximately  parallel 
to  the  surface  and  there  is  little  mixing.  In  the  turbulent 
boundary  layer  the  flow  experiences  irregular  fluctuations 
in  velocity  magnitude  and  direction  superimposed  on  the  mean 
flow,  which  remains  roughly  parallel  to  the  surface 
[Ref.  12 :p . 11] . 

The  change  from  laminar  to  turbulent  flow  is  referred  to 
as  transition.  Transition  may  be  caused  by  a  number  of 
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factors  including  disturbances  that  occur  in  the  freestream, 
surface  roughness,  changes  in  the  pressure  distribution  on 
the  surface  of  the  body,  or  an  increase  in  Reynolds  number. 
Reynolds  number  is  a  ratio  of  the  inertia  to  viscous  forces 
in  a  fluid  and  is  given  by  the  following  equation: 

Re  =  (12) 

V- 

where  U  is  the  velocity,  p  is  the  density,  1  is  the  body 
dimension,  and  p  is  the  dynamic  coefficient  of  viscosity  of 
the  fluid. 

The  presence  of  shear  forces  in  the  boundary  layer  due 
to  friction  and  the  viscous  nature  of  the  fluid  result  in 
heat  being  generated.  The  heat  is  transferred  to  the  body 
and  to  the  fluid.  The  viscous  dissipation  as  heat  in  the 
boundary  layer  results  in  an  increase  in  entropy  and 
consequently  a  decrease  in  exergy.  At  hypersonic  velocities, 
the  intense  viscous  dissipation  may  result  in  temperatures 
high  enough  to  cause  dissociation  and  ionization  of  the 
surrounding  air  [Ref.  5:p.  228].  This  phenomenon  and  its 
effects  will  not  be  discussed  in  this  paper.  However,  the 
error  introduced  by  ignoring  these  effects  can  be 
significant.  Interpolating  from  Anderson  [Ref.  5:p.  511], 
the  error,  at  M=10  and  at  30,000  meters,  in  the  ratio  T2/Tx 
across  a  normal  shock  is  approximately  24  per  cent. 
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A.  BOUNDARY  LAYER  CALCULATIONS 


An  exact  determination  of  the  flow  conditions  in  and 
surrounding  the  boundary  layer  for  a  body  of  arbitrary  shape 
is  possible  by  solving  the  three  dimensional  Navier- Stokes 
equations.  These  equations  do  not  have  a  closed  form 
solution,  but  require  elaborate  computer  methods  to  be 
solved.  However,  if  we  specify  that  the  boundary  layer 
thickness  (6)  is  much  less  than  the  length  of  the  body,  and 
that  the  Reynolds  number  is  large  (>105) ,  the  Navier- Stokes 
equations  can  be  reduced  to  the  boundary  layer  equations: 

Continuity.  HfUl  *  Hill  -  0  (13) 

ox  dy 

x  Momentum:  pu-fH+pv-fH  =  -^2*  ^  (  du}  (14) 

ox  ay  dx  dy  ay 

y  Momentum:  =0  (15) 

dy 


Energy: 


dy  dy 


u- 


dx 


dy 


(16) 


where  x  is  the  direction  parallel  to  the  surface,  y  is  the 
direction  normal  to  the  surface,  u  is  the  velocity  component 
in  the  x  direction,  v  is  the  velocity  component  in  the  y 
direction,  p  is  the  density,  pe  is  the  pressure  at  the  outer 
edge  of  the  boundary  layer,  p  is  the  viscosity  coefficient 
of  the  fluid,  h  is  the  enthalpy,  and  k  is  the  thermal 
conductivity  of  the  fluid. 
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Like  the  Navier- Stokes  equations,  the  boundary  layer 
equations  are  nonlinear.  They  are  simpler  and  easier  to 
solve,  but  still  require  computer  methods  to  obtain  exact 
solutions  in  most  cases.  Solution  methods  for  the  boundary 
layer  equations  are  presented  by  Young  [Ref.  12]  and 
Schlichting  [Ref.  14]. 

1.  Approximate  Methods 

The  purpose  of  this  thesis  was  not  to  determine  an 
exact  solution  to  the  boundary  layer  equations,  but  to  show 
that  an  exergy  analysis  can  be  conducted  on  the  boundary 
layer  and  to  provide  the  foundation  for  an  exergy  based 
design  tool  and  or  methodology.  This  concept  can  be 
demonstrated  by  making  several  approximations  that  will 
provide  rough  estimates  of  the  required  variables  without 
the  extensive  numerical  and  computer  calculations  required 
to  obtain  a  more  exact  solution  of  the  boundary  layer 
equations.  The  necessary  calculations  are  developed  in  the 
following  paragraphs. 

a.  Stagnation  Point  Calculation 

As  stated  in  the  preceding  chapter,  the  flow 
field  surrounding  a  blunt  body  traveling  at  supersonic 
velocity  will  consist  of  both  subsonic  and  supersonic  flow 
divided  by  the  sonic  line  (see  Figure  1) .  For  a  body  at  zero 
angle  of  attack,  the  flow  precisely  along  the  body  axis  will 
pass  through  the  shock  wave,  which  at  that  point  is  normal 
or  very  nearly  normal  to  the  freestream  direction,  slow,  and 
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3 top  at  the  stagnation  point.  On  a  blunt  body,  the  boundary 
layer  has  a  finite  thickness  at  the  stagnation  point.  The 
boundary  layer  thickness  at  the  stagnation  point  of  a  two 
dimensional  body  is  given  by 


6  (0) 


0.075U 

f>V 


(17) 


where  6(0)  is  the  boundary  layer  thickness  at  the  stagnation 
point,  and  V  is  the  magnitude  of  the  velocity  gradient  at 
the  stagnation  point.  [Ref.  15:pp.  208-209] 

The  velocity  gradient  is  determined  by  the 
method  of  Anderson  [Ref.  5:pp.  255-256].  It  can  be  shown 
using  Euler's  equation  applied  at  the  boundary  layer  that, 

dx  p.u.  dx 


where  the  subscript  e  indicates  conditions  at  the  edge  of 
the  boundary  layer,  u  is  the  velocity  component  in  the  x 
direction  which  is  parallel  to  the  body,  and  p  is  the 
pressure.  Assuming  a  modified  newtonian  pressure 
distribution  over  the  surface,  the  pressure  coefficient  is 
given  by  [Ref.  5:p.  53] 

Cp  =  C^sin2©,  (19) 

where  $x  is  the  angle  between  the  tangent  to  the  surface  and 
the  freestream  direction,  and  is  given  in  terms  of  y 

and  the  freestream  Mach  number  (Mx)  as 
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2 


(7  +  1  )2Ml 


(20) 


(  [ 


7  m\  47ftfJ-2(7-l) 


7  +  1 


Defining  $  as  the  angle  between  the  normal  to  the  surface 
and  the  freestream,  then  equation  19  can  be  written  as 


c°  =  c^cos2<t> 


(21) 


Defining  Cp  as 


r  _  P*-P i 

P  rr 


(22) 


where  the  subscript  1  denotes  the  freestream  condition,  then 
equation  21  can  be  written  as: 


Pa  —  =  C  cos2<|> 

ql 


(23) 


where  qx  is  the  dynamic  pressure  of  the  freestream.  Equation 
23  can  be  written  as 


Pe  9 


(24) 


Differentiating  equation  24  with  respect  to  x,  we  obtain 


Substituting  equation  25  into  equation  18,  we  have 


(25) 


du 


e  _ 


dx 


— ft** —  cos<|>sin<t> 
Peue  dx 


(26) 


Consider  the  stagnation  region,  as  shown  in  Figure  7. 
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In  this  region,  let  Ax  be  a  small  increment  of  surface 
distance  above  the  stagnation  point,  corresponding  to  small 
changes  in  <f>  of  A$.  Since  the  flow  velocity  is  so  low  in  the 
stagnation  region,  we  can  assume  [Ref.  5:p.  252] 

u*  =  (-$2>  Ax  (27) 

*  dx  s 

In  the  stagnation  region  the  following  approximations  can  be 
made  if  <f>  is  small: 


cos<|>  =  1 

(28a) 

sin<J>  =4>=A4>  =  -^ 

fv 

(28b) 

n 

Sl-S 

(28c) 

where  R  is  the  radius  of  the  nose  of  the  body  at  the 
stagnation  point. 
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(31) 

This  value  for  the  velocity  gradient  is  then  substituted 
into  the  value  for  V  in  equation  17. 

The  velocity,  temperature,  density,  Mach  number, 
and  pressure  immediately  behind  the  shock  wave  can  be 
determined  from  the  normal  and  oblique  shock  relations 
contained  in  Anderson  [Ref.  2 :pp. 67- 68 , 106] .  The  density 
behind  the  shock  along  the  stagnation  streamline  is  used  in 
equation  17  for  the  value  of  p. 

Equation  17  gives  the  boundary  layer  thickness 
at  the  stagnation  point  of  a  two  dimensional  body.  The 
boundary  layer  thickness  for  an  axisymmetric  body  will  be 
thinner  than  for  the  two  dimensional  case  [Ref.  5:p.  254] . 
While  the  difference  between  the  axisymmetric  and  two 
dimensional  cases  is  noted,  equation  17  does  provide  an 
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approximation  of  the  boundary  layer  thickness  for  the 
axisymmetric  case  at  the  stagnation  point. 

b.  Boundary  Layer  Edge  Conditions 

The  pressure  distribution  at  the  surface  of  the 
body  is  approximated  by  the  modified  newtonian  method  [Ref. 
5:pp.  53-56].  Although  the  results  provided  by  this  method 
are  not  very  accurate  at  low  supersonic  Mach  numbers,  the 
results  improve  quickly  as  Mach  number  increases  [Ref.5:pp. 
45-56] .  The  pressure  coefficient  at  the  surface  of  the  body 
using  the  modified  newtonian  method  is  given  by  equation  19. 
Recalling  from  equation  22  the  definition  for  Cp<  the 
pressure  at  the  boundary  layer  edge  over  the  nose  of  the 
body  is  approximated  by 

pe  =  CPamtq1cos20 +px  (32) 

Figure  8  shows  the  surface  pressure  distribution  over  the 
nose  of  the  hypothetical  vehicle  used  in  this  analysis  for 
selected  Mach  numbers.  The  nose  radius  (R)  in  Figure  8  is 
1.5  ft  (.4572  m)  and  the  pressure  distribution  is  shown  from 
the  stagnation  point  through  an  angle  of  68.2  degrees,  where 
the  rounded  nose  intersects  the  cone  of  the  body. 

The  total  temperature  behind  the  shock  wave 
remains  constant  everywhere  behind  the  shock.  Also,  we 
assume  that  the  total  pressure  is  constant  along  the 
streamline  that  enters  through  the  nearly  normal  portion  of 
the  shock  and  flows  along  the  boundary  layer  edge.  These 
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values  can  be  calculated  and  used  as  a  reference  to  obtain 


local  values  of  temperature  and  pressure. 


The  total  pressure  in  the  freestream  is  given  by 
P0  =  P[(i  +  2lJ:Mf)  V3!]  (33) 

where  P01  is  the  stagnation  pressure  in  the  freestream,  P  is 
the  static  pressure  of  the  freestream.  For  these  analyzes 
all  freestream  atmospheric  conditions  were  taken  from 
Reference  [6]  at  an  altitude  of  30,000  meters.  The  ratio  of 
stagnation  pressures  across  the  shock  is 
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(34) 


e 


s ,  -  s, 


R' 


where  R'  is  the  specific  gas  constant,  and  Sj-Sj^  is  given  by 


s2-s1  =  CplnKl  +  il  - 1) ) 


2*  (7-1  )Mf 


-Rln[l  +  JiL  (Ml  -1)  ] 


)  3 


(35) 


where  cp  is  the  specific  heat  at  constant  pressure  of  the 
gas.  The  stagnation  pressure  behind  the  shock  is  then 
obtained  by  taking  the  product  of  equations  33  and  34 


(36) 


The  stagnation  temperature  in  the  freestream  is 
calculated  by 

T0i  =  T{l  +  (37) 

where  T  is  the  freestream  temperature.  The  total  temperature 
across  a  normal  shock  is  constant,  therefore 


T0i  =  T0j  (38) 

If  we  substitute  T02  for  T01  in  equation  37,  and  with  a 
known  Mach  number,  the  local  temperature  can  be  determined. 

The  sonic  point  on  the  body  nose  is  located  with 
the  same  method  used  in  determining  the  shock  shape,  and  is 
described  in  Appendix  A.  Since  the  Mach  number  at  the  sonic 
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line  is  one,  Che  local  pressure  at  the  sonic  point  can  be 
obtained  by  substituting  the  value  of  P02,  from  equation  36, 
for  P01  in  equation  33.  Once  the  sonic  point  is  located,  the 
angle  between  the  body  centerline  and  the  sonic  point  can  be 
easily  calculated.  This  angle  will  be  designated  the  sonic 
point  angle  (SPA)  for  calculation  purposes.  The  total  angle 
that  the  flow  makes  as  it  rounds  the  nose  is  68.2  degrees. 
Subtracting  the  SPA  from  68.2  gives  the  remaining  angle, 
designated  A1#  that  the  flow  must  turn  through  in  traveling 
from  the  sonic  line  to  the  cone  portion  of  the  body. 

Treating  angle  Ax  as  a  Prandtl -Meyer  expansion  angle,  the 
pressure,  Mach  number,  and  temperature  of  the  flow 
downstream  of  the  expansion  can  be  calculated. 

Before  proceeding  further,  it  is  necessary  to 
introduce  the  Prandtl-Meyer  function  [Ref.  2:p.  134],  which 
is  given  the  symbol  v 


v  (M) 


111  tan"1  lZl(M2- 1)  -tan"VM2-l 
7-1  y  +1 


(39) 


As  can  be  seen  from  equation  39,  the  Prandtl-Meyer  function 
is  a  function  of  the  Mach  number  and  y.  Let  v  (M-^)  be  the 
Prandtl-Meyer  function  before  the  expansion,  and  MM2)  be 
the  Prandtl-Meyer  function  after  the  expansion.  Then  the 
angle  A1  that  the  flow  turns  through  is 
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At  the  sonic  line  the  Mach  number  is  one,  and 
therefore  i»(M1>  in  equation  40  is  zero  and  Ar  MM2)  .  Putting 
the  known  value  of  Ax  into  equation  39  for  v  ( M) ,  the  Mach 
number  after  the  expansion  can  be  determined.  This  process 
is  repeated  at  the  point  wb"?re  the  cone  joins  the  body.  In 
this  fashion  the  flow  Mach  number  along  the  length  of  the 
cylindrical  portion  of  the  vehicle  can  be  found,  with  the 
Mach  number  known  and  the  previously  calculated  values  of 
total  pressure  and  total  temperature,  the  local  pressure  and 
temperature  along  the  cylindrical  portion  of  the  body  can  be 
determined  using  equations  33  and  37.  These  values  will  be 
taken  as  the  conditions  at  the  edge  of  the  boundary  layer, 
c.  Boundary  Layer  Thickness 

The  boundary  layer  is  formed  as  a  result  of  the 
no  slip  condition  that  exists  at  the  wall  of  a  body  immersed 
in  a  flow  (see  Figure  6) .  As  the  flow  progresses  along  the 
length  of  the  body,  the  thickness  (5)  of  the  layer  of 
retarded  flow  increases  as  more  and  more  fluid  becomes 
affected  [Ref.  14 :p.  25] .  The  thickness  of  the  boundary 
layer  on  a  flat  plate  at  zero  incidence  can  be  calculated 
using  the  method  developed  by  Blasius  [Ref.  16].  For  laminar 
flow  the  boundary  layer  thickness  was  found  to  be 


where  x  is  the  distance  along  the  plate,  and  Rex  is  the 
Reynolds  number  (see  equation  12)  as  a  function  of  distance 
along  the  plate.  For  a  turbulent  boundary  layer  the 
thickness  is  given  by  [Ref.  l2:p.  13] 


0.37x 
(Rex) l/s 


(42) 


As  previously  stated,  the  boundary  layer  on  an 
axisymmetric  body  is  thinner  than  on  a  flat  plate.  However, 
as  a  first  approximation,  the  boundary  layer  over  a  flat 
plate  with  a  length  equal  to  the  distance  traveled  by  the 
flow  from  the  stagnation  point  to  the  tail  end  of  the  body 
will  be  used  to  approximate  the  boundary  layer  thickness. 
This  approximation  is  reasonable  when  6<<R  as  it  is  in  this 
case. 

Since  there  is  a  dramatic  difference  in  the 
growth  of  laminar  vs  turbulent  boundary  layers,  it  is 
necessary  to  establish  the  point  where  transition  occurs.  In 
actuality  transition  may  occur  over  a  finite  distance  making 
prediction  of  a  point  of  transition  extremely  difficult. 
Research  to  determine  exactly  how  transition  occurs  and 
methods  to  predict  its  occurrence  are  ongoing.  A  formula 
developed  by  Anderson  [Ref.  5:p  280]  to  predict  the  point  of 
transition  will  be  used  in  this  analysis. 
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(43) 


log10  (ReT)  =  6.421  exp  [1. 209x10  ~4 Mg  • 64 1 ) 

where  ReT  is  the  Reynolds  number  at  which  transition  occurs, 
and  Me  is  the  Mach  number  at  the  boundary  layer  edge.  It 
should  be  recognized  that  the  data  used  to  develop  this 
equation  was  unknown  to  the  author  and  that  equation  42  was 
being  used  as  an  approximation. 

B.  BASE  PRESSURE  DETERMINATION 

If  the  assumption  is  made  that  the  launch  vehicle  is  a 
projectile  with  no  exhaust  plume,  the  pressure  behind  the 
body  in  the  wake  is  assumed  to  be  the  freestream  static 
pressure.  This  assumption  is  supported  by  Anderson  [Ref. 

5 : pp  117-138]  who  shows  that  for  a  cylindrical  body  with 
fineness  ratio  greater  than  six,  the  pressure  asymptotically 
approaches  freestream  pressure.  The  assumption  of  no  exhaust 
plume  allows  calculation  of  the  entropy  change  as  a  result 
of  the  boundary  layer  on  the  body,  but  ignores  the  effects 
of  boundary  layer/plume  interaction. 

With  the  pressure  in  the  wake  being  approximately  equal 
to  static  pressure  and  the  stagnation  pressure  behind  the 
shock,  from  equation  36,  the  Mach  number  in  the  wake  can  be 
determined  using  equation  33.  The  difference  between  the 
Mach  number  in  the  wake  and  in  the  flowfield  can  be  used  to 
calculate  the  difference  between  the  velocities  (U-u)  needed 
to  solve  equation  6.  These  values  are  summarized  in  Table  3. 
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C.  WAKE  PROFILE 


The  wake  boundary  is  established  from  the  body  plus  the 
boundary  layer  thickness  at  the  base  of  the  body.  Traveling 
downstream  for  several  body  diameters  and  assuming  the  same 
rate  of  growth  of  the  boundary  layer,  the  distance  from  the 
centerline  to  the  point  where  the  wake  velocity  approaches 
the  freestream  velocity  can  be  determined.  This  is  shown  in 
Figure  9.  For  calculation  purposes,  a  distance  of  five 
diameters  downstream  of  the  base  was  arbitrarily  chosen. 


TABLE  3:  COMPARISON  OF  FREESTREAM  AND  WAKE  MACH  NUMBER 


Freestream  Mach  (M^) 

Wake  Mach(M„) 

1.1 

1.0991 

.  ,  3  92 

2.0 

1.7879 

.8939 

3.0 

2.2769 

.7590 

4.0 

2.6351 

.6588 

5.0 

2.9218 

.5844 

6.0 

3.1630 

.5272 

7.0 

3.3725 

.4818 

o 

• 

CD 

3.5586 

.4448 

9.0 

3.7265 

.4141 

10.0 

3.8800 

.3880 

Figure  9:  Illustration  of  Wake  Profile. 
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The  variation  of  u  in  the  wake  from  the  centerline  to 
the  point  where  u»U  is  given  by 

u  =  Uu  *  Ux  (44) 

where  Uw  is  the  minimum  wake  velocity,  Ux  is  the  freestream 
velocity,  and  x  is  a  scaling  variable  where  Osxsl.  The 
distance  from  the  centerline  to  the  point  where  u-U  is  given 
by 

y  =  (**  +  8C)  n  (45) 

where  Rb  is  the  radius  of  the  base  of  the  body  and  <5t  is  the 
thickness  of  the  boundary  layer  extended  to  the  point  in  the 
wake  where  the  profile  is  being  calculated. 

D .  DRAG  DETERMINATION 

The  drag  of  the  body  is  calculated  using  equation  6 
which  is  repeated  here: 

D  =  f  pu(U-u)dA  (46) 

J  A 

Since  the  velocity  in  the  wake  varies  along  the  centerline, 
designated  the  x  axis,  and  with  the  distance  from  the 
centerline,  designated  the  y  axis,  equation  46  becomes  the 
following  double  integral: 

D  -  fyf12nrpu(U-u)dr  (47) 

where  u  and  y  are  given  by  equations  44,  and  45 
respectively.  The  density  in  equation  47  is  the  density  at 
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the  location  in  the  wake  where  the  profile  is  determined. 

For  these  calculations  this  value  was  approximated  using  the 
perfect  gas  equation.  The  pressure  is  taken  as  freestream 
and  the  temperature  is  that  calculated  for  the  boundary- 
layer  temperature  just  forward  of  the  base  of  the  body 
ignoring  plume  effects. 

E.  ENTROPY  CALCULATION 

With  the  drag  determined,  the  entropy  flux  increase  as  a 
result  of  the  boundary  layer  can  now  be  calculated  using 
equation  11  which  is  repeated  here 

DU  =  71  — —  (48) 

1  sec 

where  U  and  Tx  are  the  freestream  velocity  and  temperature 
respectively,  and  D  is  the  drag. 

For  the  body  shape  used  in  this  analysis  and  with  the 
approximations  used  to  calculate  the  required  flowfield 
parameters,  the  entropy  increase  is  presented  in  Table  4. 

As  can  be  seen  from  Table  4,  the  entropy  increase  in  the 
boundary  layer  starts  to  decrease  above  M=7.  This  occurs  due 
to  the  boundary  layer  becoming  thinner  at  higher  Ma  li 
numbers.  The  thinner  boundary  layer  is  a  result  of  the 
transition  point  moving  aft  on  the  body.  Initially  the 
transition  point  moves  forward  toward  the  nose  of  the  body. 
Between  6sMs7,  the  transition  point  stops  moving  toward  the 
nose  and  starts  moving  back  toward  the  tail.  This  is 
occurring  as  a  result  of  the  laminar  to  turbulent 
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TABLE  4:  ENTROPY  INCREASE  DUE  TO  THE  BOUNDARY  LAYER 


Mach  number 

Rate  of  Entropy  increase 
(KJ/sec) 

1.1 

2918 

2 

21230 

3 

48316 

4 

70635 

5 

84925 

6 

92074 

7 

92671 

8 

87442 

9 

77123 

10 

62372 

transition  point  as  calculated  by  equation  43.  The 
transition  Reynolds  number  increases  exponentially  with 
increasing  Mach  number.  For  a  given  transition  Reynolds 
number,  the  point  of  transition  is  obtained  by  substituting 
ReT  for  Re  and  x  for  1  in  equation  12,  and  rearranging 


x 


Re 

~PV~ 


(49) 


As  ReT  increases  exponentially,  the  numerator  in  equation  49 
grows  faster  than  the  denominator,  as  shown  in  Figure  10. 

The  result  of  this  is  that  the  transition  point  moves  toward 
the  tail  of  the  body  and  the  boundary  layer  remains  laminar 
longer.  Since  entropy  is  produced  at  a  higher  rate  in 
turbulent  boundary  layers,  the  net  result  is  that  less 
entropy  production  occurs. 
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Mach  Number. 


This  phenomenon  is  a  result  of  using  equation  43  to 
predict  the  transition  point  and  may  not  be  an  accurate 
representation  of  actual  flow.  As  previously  stated,  the 
data  that  equation  43  was  developed  from  was  unknown  to  the 
author,  and  the  conditions  of  this  study  may  be  outside  the 
range  for  which  equation  43  is  valid. 
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IV.  CONCLUSIONS 


Recall  from  the  exergy  equation  (equation  1)  that  an 
increase  in  entropy  results  in  a  corresponding  decrease  in 
exergy.  A  primary  advantage  of  exergy  analysis  is  the 
identification  of  those  processes  or  areas  that  have  the 
highest  exergy  loss  and  hence  the  greatest  potential  for 
improved  efficiency.  A  method  to  quantify  the  exergy 
decrease  resulting  from  entropy  production  across  the  shock 
wave  and  in  the  boundary  layer,  ignoring  high  temperature 
effects,  has  been  presented. 

It  has  been  demonstrated  that  the  exergy  decrease  due  to 
the  production  of  entropy  across  the  shock  wave  and  in  the 
boundary  layer  can  be  calculated  for  a  typical  space  launch 
vehicle.  This  information  can  be  used  to  highlight  those 
areas  which  have  the  greatest  potential  for  improved 
efficiency. 

For  the  conditions  of  this  study  and  within  the  accuracy 
of  the  approximations,  it  is  clearly  evident  that  entropy 
production  is  much  greater  in  the  boundary  layer  than  across 
the  shock  wave  for  a  given  Mach  number.  This  indicates  that, 
all  other  factors  being  equal,  design  efforts  should 
concentrate  on  reducing  the  entropy  production  in  the 
boundary  layer. 
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V.  RECOMMENDATIONS  POR  FURTHER  RESEARCH 


This  thesis  applies  the  method  of  exergy  analysis  to  two 
areas  of  interest  for  space  launch  vehicles.  In  doing  so,  it 
provides  an  initial  step  in  the  development  of  an  exergy 
methodology  tool  for  design  purposes. 

Following  the  work  of  Czysz  [Ref.  21] ,  application  of 
this  methodology  to  launch  vehicle  propulsion  systems  would 
be  beneficial. 

Alternately,  variations  in  launch  vehicle  design  could 
be  examined  with  the  goal  of  finding  the  optimum  design  to 
minimize  entropy  production  across  the  shock  wave  and  within 
the  boundary  layer. 

It  would  also  be  useful  to  determine  the  magnitude  of 
the  error  for  the  approximations  used  in  this  study.  This 
could  be  accomplished  by  comparing  the  results  with  those 
obtained  from  exact  flowfield  calculations  using  CFD 
methods.  Also  the  effects  of  the  boundary  layer/plume 
interactions  on  entropy  production  could  be  investigated. 
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APPENDIX  A 
THE  SHOCK  WAVE  MODEL 


This  program  computes  the  entropy  increase  across  the 
shock  wave  of  a  cylindrical  space  launch  vehicle  with  a 
blunt  cone  nose.  The  program  is  written  for  Mach  numbers  of 
1.1  to  10,  and  for  atmospheric  conditions  at  an  altitude  of 
30,000  meters.  This  program  is  written  in  Fortran  77  code. 
Comment  statements  are  preceded  by  a  "C" . 


C  Variable  declarations 

M  =  Freestream  Mach  number 
BETA  »  Cotangent  of  the  Mach  number 

THETA  =  The  shock  wave  angle  for  sonic  flow  behind  the 
shock 

DELTA  -  The  maximum  stream  deflection  angle  that  can 
occur  without  shock  separation 
NU  =  Angle  between  the  sonic  line  and  a  line  normal  to 
the  freestream  direction 
GAMMA  -  Ration  of  specific  heats  (1.4  for  air) 

SIGMA  *  Angle  of  shock  wave  at  the  center  of  mass 
streamline 

A  =  Area  ratio  for  isentropic  flow 
P  =  Stagnation  pressure  ratio 

RATIO  »  Ratio  of  the  distance  from  the  x-axis  to  the 

shock  sonic  point  to  the  distance  from  the  x-axis 
to  the  body  sonic  point 

YSB  *  Distance  from  the  centerline  of  the  body  to  the 
sonic  point  on  the  body 

X0  -  Distance  from  foremost  point  of  detached  shock  to 
intercept  of  its  asymptote  on  the  x-axis 
Y  »  Coordinate  perpendicular  to  freestream  flow 
X  a  Coordinate  in  the  direction  of  the  freestream 
I , J, Z  ■  Increment  counters 
PI  a  Radian  equivalent  of  180  degrees 

DIFF  -  The  difference  between  the  shock  wave  angle  and 
the  asymptote  angle  of  the  hyperbola  that  is 
defined  by  the  shock  wave 

CANGLE  »  The  constant  angle  between  the  asymptote  of 
the  hyperbola  and  the  freestream  direction 


40 


ooonnn 


SANGLE  =  The  angle  between  the  freestream  direction  and 
the  tangent  to  the  shock  wave 
ANGLED  =  The  same  angle  as  "SANGLE"  but  in  degrees 
instead  of  radians 

TOTALS  =  The  total  change  in  entropy  across  the  shock 
wave 


C 

C 

C 

C 

C 

C 

C 

C 


START  PROGRAM 
INTEGER  I,  J,  Z 
PARAMETER  (Z=10) 

REAL  M(Z)  ,  BETA (Z) ,  THETA ( Z)  ,  DELTA ( Z) , 

+  NU (Z) ,  PI,  YSB(Z),  SIGMA ( Z)  ,  A{Z),  P(Z), 

+  RATIO ( Z) ,  X0(Z),  X (7500) ,  Y{7500),  DIFF,  B(Z), 

+  CANGLE { Z ) ,  SANGLE (7500) ,  ANGLED ( 7500 ) ,  TOTALS 

PI=3 . 14159 
M ( 1) =1.1 
M { 2 ) =2 

BETA  (1)  =SQRT  (M  (1)  **2 .0-1.0) 

BETA (2) =SQRT(M(2) **2 .0-1.0) 

DO  05  1=3, Z 

M ( I) =M( I- 1) +1 . 0 

BETA (I) =SQRT (M ( I) **2.0 -1.0) 

05  CONTINUE 

CALL  SHOCK (Z,M, THETA) 

CALL  MAX (Z,M, DELTA) 

CALL  AVG ( Z , M, NU) 

CALL  DISYSB (Z, DELTA, YSB) 

CALL  CENTER (Z, THETA, BETA, SIGMA) 

CALL  AREA(Z,M, A) 

CALL  PRESS (Z, SIGMA, M, P) 

CALL  SPRATIO ( Z, P, A, NU, RATIO) 

CALL  DISXO ( Z , YSB , BETA, RATIO , THETA, X0 ) 

The  following  do  loop,  for  each  Mach  number,  calculates 
the  angle  of  the  asymptote  of  the  hyperbola  that 
defines  the  shock.  Then  the  point  on  the  x  and  y 
axis  where  the  difference  between  the  shock  wave 
angle  and  the  asymptote  angle  is  less  than  1  per 
cent  is  found.  This  defines  the  limit  of  the  shock  wave 
for  purposes  of  calculating  the  entropy  change  across 
the  shock.  The  change  in  entropy  is  then  calculated. 

DO  15  J=l,Z 

B ( J) =X0 { J) /BETA ( J) 

CANGLE  (J)  =  ATAN  { B  ( J )  /X0  (J)  ) 

CALL  SHAPE (X0 ( J) , BETA ( J) , X, Y, ) 

CALL  SWANGLE ( X0 ( J ) , BETA ( J) , Y, SANGLE , ANGLED ) 

DO  10  1=2,7500 

IF  (SANGLE (I) .EQ.0)  GO  TO  10 
DIFF= (SANGLE (I) -CANGLE (J) ) /CANGLE (J) 

IF (DIFF. GT. 0.01)  GO  TO  10 
IF (DIFF. LE. 0.01)  THEN 

CALL  TOTAL (M(J) , X0 (J) , BETA( J) ,X(I) , TOTALS) 

GO  TO  15 
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END  IF 


10  CONTINUE 
15  CONTINUE 
STOP 
END 


C  Calculate  the  shock  wave  angle  for  sonic  flow  behind 
C  the  shock 

INTEGER  I,Z 

REAL  M ( Z ) ,  THETA ( Z ) ,  GAMMA 
GAMMA= 1 . 4 
DO  05  1=1, Z 

THETA ( I ) =ASIN (SQRT ( (1.0/4. Q*GAMMA*M ( I ) **2 . 0 ) ) 

+  *( (GAMMA+1.0) *M(I) **2.0- (3.0-  GAMMA) 

+  + (SQRT ( (GAMMA+1.0) * ( (GAMMA+1.0) *M(I) **4 . 0-2 . 0 

+  * (3 .0 -GAMMA) *M(I) **2.0+ (GAMMA+9 . 0) )))))) 

05  CONTINUE 
RETURN 
END 

C  Calculate  the  maximum  deflection  angle,  in  radians, 

C  for  a  given  freestream  Mach  number  in  air.  For  this 

C  calculation,  the  maximum  deflection  angle  without  shock 

C  detachment  was  calculated  from  NACA  report  1135,  chart 
C  5.  An  equation  for  the  line  of  maximum  deflection  was 

C  obtained  by  performing  a  curve  fit  of  the  data 

C  presented  in  the  chart. 

SUBROUTINE  MAX (Z,M, DELTA) 

C  NLOG  =  Natural  log  of  Mach  number  between  l.i  and  10 

C  TVAR  =  Natural  log  of  shock  wave  detachment  angle  in 

C  degrees 

INTEGER  I,Z 

REAL  M (Z) ,  NLOG (10),  TVAR(10),  DELTA { Z ) 

DO  05  1=1, Z 

NLOG ( I ) =  LOG (M ( I) ) 

TVAR (I) =3. 80568+0. 41631*L0G (NLOG (I) ) 

DELTA(I) = (EXP (TVAR (I) ) )/57. 29578 
05  CONTINUE 
RETURN 
END 


C  Calculate  the  angle  between  the  sonic  line  and  a  line 
C  normal  to  the  freestream  direction 
SUBROUTINE  AVG(Z,M,NU) 

INTEGER  I,Z 
REAL  M ( Z) ,  NU (Z) 

DO  05  1=1, Z 

NU(I) =(3.9901+56. 9298*LOG(M( I) ) )/57. 29578 
05  CONTINUE 
RETURN 
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END 

C  Calculate  the  distance  from  the  body  centerline  to  the 
C  sonic  point  on  the  body 

SUBROUTINE  DISYSB ( Z , DELTA, YSB ) 

C  DISTl=Distance  along  a  line  drawn  tangent  to  the  body 
C  sonic  point  from  the  x-axis  to  the  sonic  point  on 

C  the  body 

C  R=Radius  of  the  body  apex  (in  feet) 

C  ANG3=The  angle  between  the  x-axis  and  a  line 
C  perpendicular  to  DIST1  at  the  body  sonic  point  and 

C  extending  from  the  sonic  point  to  the  x-axis 

INTEGER  I,Z 

REAL  PI,  YSB(Z),  DISTl(lO),  R,  ANG3(10),  DELTA (Z) 

PI=3 .14159 
R-1.5 

DO  05  1=1, Z 

ANG3 (I) =PI- (PI/2.0) - DELTA ( I ) 

DIST1 ( I ) =R* ( (SIN (ANG (I) ) /SIN (DELTA (I) ) ) ) 

YSB ( I) =DIST1 ( I) *  S IN ( DELTA ( I ) ) 

05  CONTINUE 
RETURN 
END 

C  Calculate  the  angle  of  the  shock  wave  at  the  center 
C  of  mass  streamline  between  the  shock  wave  vertex 

C  and  the  sonic  point  of  the  shock  wave 

SUBROUTINE  CENTER ( Z , THETA, BETA, SIGMA) 

INTEGER  I,  Z 

REAL  THETA (Z) ,  BETA(Z) ,  SIGMA(Z) 

DO  05  1=1, Z 

SIGMA (I) =ATAN ( ( ( TAN ( THETA ( I ) ) ) **2 . 0* (3 . 0/2 . 0) **2 . 0 
+  - (1.0/BETA (I) ) **2.0) * ( (3 .0/2.0) **2.0 -1.0) ) **0.5) 

05  CONTINUE 
RETURN 
END 

C  Calculate  the  area  ratio  for  isentropic  flow 
SUBROUTINE  AREA(Z,M, A) 

INTEGER  I,Z 

REAL  GAMMA,  M(Z),  A(Z) 

GAMMA=1 . 4 
DO  05  1=1, Z 

A (I) =( (GAMMA+1 .0) /2 .0) ** ( (GAMMA+1.0) / (2.0 
+  *GAMMA- 1.0) ) ) *M(I) * (1.0+ ( (GAMMA- 1 . 0 ) /2 . 0 ) 

+  *M(I) **2.0) ** (- ( (GAMMA+1.0) / (2 .0* (GAMMA-1.0) ) ) ) 

05  CONTINUE 
RETURN 
END 
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C  Calculate  the  ratio  between  stagnation  pressure  before 

C  and  after  the  shock  wave  for  a  given  shock  wave  angle 

SUBROUTINE  PRESS ( Z , SIGMA, M  ?) 

INTEGER  I ,  Z 

REAL  GAMMA,  SIGMA(Z),  M(Z),  P(Z) 

GAMMA  =  1.4 
DO  05  1=1, Z 

PI  (I) =1.0/ {  (  ( (GAMMA+1 .Q)*M(I)**2.0* 

+  SIN(SIGMA(I) ) **2.0) /( (GAMMA-1.0) * 

+  M ( I ) **2 . Q*SIN (SIGMA ( I ) ) **2 . 0+2 . 0) ) ** (GAMMA/ 

+  (GAMMA- 1.0) )*( ( GAMMA+1. 0)/ (2. 0*GAMMA*M(I)* *2.0* 

+  SIN(SIGMA(I) ) **2.0- (GAMMA-1.0) ) )** 

+  (1.0/ (GAMMA-1.0) ) ) 

05  CONTINUE 
RETURN 
END 

C  Calculate  the  ratio  of  the  distance  from  the  x-axis  to 
C  the  shock  sonic  point  to  the  distance  from  the  x- 
C  axis  to  the  body  sonic  point 

SUBROUTINE  SPRATIO ( Z , P , A, NU, RATIO) 

INTEGER  I,Z 

REAL  P(Z),  A ( Z j  ,  NU(Z),  RATIO 
DO  05  1=1, Z 

RATIO ( I) =  1.0/ ( ( 1 . 0 - P (I ) *A (I) *COS (NU ( I ) ) )**0.5> 

05  CONTINUE 
RETURN 
END 

C  Calculate  the  distance  from  the  foremost  point  of  the 

C  detached  shock  wave  to  the  intercept  of  its  asymptote 

C  on  the  x-axis 

SUBROUTINE  DISX0 ( Z , YSB , BETA, RATIO, THETA, X0 ) 

INTEGER  I ,  Z 

REAL  YSB(Z),  BETA ( Z ) ,  RATIO(Z),  THETA ( Z) ,  X0(Z) 

DO  05  1=1, Z 

X0 (I) =YSB (I) *BETA(I) *RATIO (I) * ( ( BETA ( I ) )**2.0* 

+  (TAN (THETA (I) ))**2.0-1.0)**0.5 

05  CONTINUE 
RETURN 
END 

C  Calculate  the  x  and  y  coordinates  for  any  point  on  the 
C  shock  wave 

SUBROUTINE  SHAPE (X0 , BETA, X, Y) 

INTEGER  I 

REAL  X (7500) ,  Y(7500),  BETA,  X0 
X(l) =0.1 
Y (1) =0 

DO  05  1=2,7500 

X(I) =X(I-1) +0.2 
IF (X (I) . LT.X0)  THEN 


Y(I) =0 

ELSEIF  (X(I) .GE.XO)  THEN 
Y{I) =SQRT (X (I ) **2.0-X0**2 .0) /BETA 
END  IF 
05  CONTINUE 
RETURN 
END 

C  Calculate  the  angle  between  the  freestream  and  a  line 

C  drawn  tangent  to  the  shock  wave  at  any  point 

SUBROUTINE  SWANGLE (X0 , BETA, Y, SANGLE, ANGLED) 

INTEGER  I 

REAL  X0,  BETA,  Y(7500),  SANGLE { 7500 ) ,  PI,  ANGLED ( 7500 ) 
PI-3 .14159 
DO  05  1=1,7500 

IF  (Yd).EQ.O)  THEN 
SANGLE (I) =0 

ELSEIF  (Y(I).GT.O)  THEN 

SANGLE (I )=ATAN( (SQRT (X0**2 . 0+BETA**2 . C*Y d) **2.0)  ) 
+  / (BETA**2.0*Y{I)  )  ) 

ANGLED (I) =SANGLE (I) *180/PI 
END  IF 
05  CONTINUE 
RETURN 
END 

Calculate  the  total  change  in  entropy  across  the  shock 
wave.  All  calculations  are  in  SI  units 
RO=The  density  of  air.  For  these  calculations  this 

value  is  taken  from  the  ICAO  standard  atmosphere 
tables  at  an  altitude  of  30,000  meters. 

TEMP=The  temperature  in  degrees  Kelvin  of  the 

atmosphere  at  30,000  meters  altitude,  taken  from 
the  ICAO  standard  atmosphere  tables 
GAMMA=Ratio  of  specific  heats 

SGC=Specific  gas  constant.  Assuming  an  ideal  gas,  for 
air  this  value  is  287  Joules/Kg-K 
SSA=Local  speed  of  sound  in  air 

VEL=Velocity,  in  meters  per  second,  of  the  freestream 
SHCV=Specif ic  heat  of  air  at  constant  volume 
PI=Radian  equivalent  of  180  degrees 

LLINT=Lower  limit  of  integration.  The  starting  point, 
located  at  the  shock  wave  vertex,  for  numerical 
integration  to  calculate  the  entropy  change  across 
the  shock  wave 

ULINT=Upper  limit  of  integration.  The  ending  point  for 
numerical  integration  across  the  shock  wave 
H=The  integration  step  size  for  application  cf 
Simpson's  rule 

N=The  number  of  Simpson's  rule  integrations,  and  the 
ending  point  for  summation  of  individual 
integrations 
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CENTERX=The  center  point,  along  the  x-axis,  for  each 
Simpson's  rule  integration 
TSTEP-Specif ies  the  beginning  (TSTEP(l)},  center 
(TSTEP (2) ) ,  and  end  (TSTEP<3))  points  for 
Simpson's  rule  application  for  each  individual 
summation 

RADIUS=The  distance  from  the  centerline  x-axis  to  the 
shock  wave  at  each  TSTEP  location 
TANGLE =The  shock  wave  angle  at  each  point  of 

intersection  of  the  radius  and  the  shock  wave 
DELTAE=Change  in  entropy  per  mass  flow  across  the  shock 
at  each  point  of  intersection  of  RADIUS  with  the 
shock  wave 

DELTAS=Distance  along  the  shock  at  each  interval  of 
integration 

HOLD=Temporary  storage  variable 
INTEGER  I,  K 

REAL  XO,  X,  BETA,  DELTAE ( 3 ) ,  TOTALS,  LLINT,  ULINT,  H, 

+  CENTERX ,  PI,  RADIUS ( 3 ) ,  SHCV,  GAMMA,  TANGLE ( 3 ) , 

+  TSTEP (3),  RO,  VEL,  SSA,  TEMP,  M,  DELTAS 

R0=. 01841 
TEMP=226 . 509 
GAMMA= 1 . 4 
SGC=287 

SSA-SQRT  ( GAMMA*  SGC*  TEMP) 

VEL=M*SSA 
SHCV=717.5 
PI  =  3. 14159 
LLINT-XO* .3048 
ULINT=X* .3048 
HOLD=0 
H=0 . 5 

N=INT ( (ULINT-LLINT) / (2.0*H) ) 

CENTERX=LLINT+ . 05 
DO  05  I-LLINT,  N 
DO  10  K-1,3 
TSTEP ( 1 ) =CENTERX -.05 
TSTEP (2) -CENTERX 
TSTEP ( 3 ) = CENTERX + . 05 
IF  {TSTEP { K) .LE. LLINT)  THEN 
RADIUS (K)=0 
TANGLE (K) =0 
DELTAE (K) =0 
GO  TO  10 
END  IF 

RADIUS {K) - (SQRT (TSTEP (K) **2 . 0- LLINT* *2 . 0) ) /BETA 
TANGLE (K) =ATAN ( (SQRT (LLINT* *2 . 0+BETA**2 . 0*RADIUS (K) 

+  **2.0) ) / ( BETA* *2 . 0*RADIUS (K) ) ) 

DELTAE (K)=SHCV* (LOG ( (2 . 0*GAMMA*M**2 . 0*SIN (TANGLE (K) ) 

+  **2 .0- (GAMMA-1.0)  )./ (GAMMA+1.0)  ) -GAMMA* 

+  LOG ( ( (GAMMA+1.0) *M**2 . 0*SIN (TANGLE (K) ) **2.0) / 

+  ( (GAMMA- 1.0) *M**2 .0*SIN (TANGLE (K) ) **2. 0+2.0) ) ) 
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10  CONTINUE 

DELTAS  ==  (0.1*.  3048)  /COS  (TANGLE  (2)  ) 

HOLD=HOLD+ ( ( (H/3.0) * (DELTAE ( 1 ) +4 . 0*DELTAE (2) 

+  +DELTAE ( 3 ) ) ) *RO*2 . 0*PI*RADIUS (2) *DELTAS*VEL* 

+  SIN (TANGLE (2) ) ) 

CENTERX=CENTERX+0 . 1 
05  CONTINUE 

TOTALS =HOLD*TEMP 

RETURN 

END 


47 


LIST  OF  REFERENCES 


1.  Ahern,  J.E.,  The  Exergy  Method  of  Energy  Systems 
Analysis,  John  Wiley  and  Sons,  1980. 

2.  Anderson,  J.D.,  Modern  Compressible  Flow  With  Historical 
Perspective,  McGraw  Hill  Book  Co.,  1990. 

3.  National  Advisory  Committee  for  Aeronautics  Technical 
Note  1921,  Approximate  Method  for  Predicting  Form  and 
Location  of  Detached  Shock  Waves  Ahead  of  Plane  of 
Axially  Symmetric  Bodies,  by  W.E.  Moeckel,  1949. 

4.  National  Advisory  Committee  for  Aeronautics  Technical 
Note  4170,  A  Reexamination  of  the  Use  of  Simple  Concepts 
for  Predicting  the  Shape  and  Location  of  Detached  Shock 
Waves,  by  E.S.  Love,  December  1957. 

5.  Anderson,  J.D.,  Hypersonic  and  High  Temperature  Gas 
Dynamics,  McGraw  Hill  Book  Co.,  1989. 

6.  International  Civil  Aviation  Organization,  Standard 
Atmosphere  Manual,  Government  Printing  Office, 
Washington,  DC,  1964. 

7.  National  Advisory  Committee  for  Aeronautics  Report  1135, 
Equations,  tables,  and  Charts  for  Compressible  Flow,  by 
Ames  Research  Staff,  1953. 

8.  Sears,  W.R.,  General  Theory  of  High  Speed  Aerodynamics, 
Princeton  University  Press,  1954. 

9.  Crocco,  L.,  Eine  neue  Stromfunktion  fur  die  der  Bewegung 
der  Gase  mi t  Rotation,  Z.  Angew.  Math.  Mech.  vol.  17, 
1937,  pp.  1-7. 

10.  Oswatitsch,  K.,  Der  Luftwi der stand  als  Integral  des 
Entropiestromes ,  Nachr.  Ges.  Wiss.  Gottingen.  Math. 

Phys.  Kl.,  1945. 

11.  Prandtl,  L.,  Uber  Flussigkeitsbewegung  bei  sehr  Kleiner 
Reibung,  Proc  III  International  Mathimatical  Congress, 
Heidelberg,  Germany,  1904. 

12.  Young,  A.D.,  Boundary  Layer,  American  Institute  of 
Aeronautics  and  Astronautics,  1989. 

13.  Anderson,  J.D.,  Fundamentals  of  Aerodynamics,  McGraw 
Hill  Book  Co.,  1984. 


48 


14.  Schlichting,  H. ,  Boundary  Layer  Theory,  7th  Ed.,  McGraw 
Hill  Book  Co. ,  1979 . 

15.  Moran,  J.,  An  Introduction  to  Theoretical  and 
Computational  Aerodynamics,  John  Wiley  and  Sons,  Inc., 
1984. 

16.  National  Advisory  Committee  for  Aeronautics  Technical 
Memorandum  1256,  The  Boundary  Layers  in  Fluids  With 
Little  Friction,  by  H.  Blasius,  1950. 

17.  von  Karman,  T.,  and  Moore,  N.B.,  Resistance  of  Slender 
Bodies  Moving  With  Supersonic  Velocities  With  Special 
Reference  to  Projectiles,  A.S.M.E.,  Transactions,  pp. 
303-310,  1932. 

18.  National  Advisory  Committee  for  Aeronautics  Technical 
Note  2137,  An  Analysis  of  Base  Pressure  at  Supersonic 
Velocities  and  Comparison  With  Experiment,  by  D.R. 
Chapman,  July  1950. 

19.  National  Advisory  Committee  for  Aeronautics  Technical 
Note  3393,  An  Experimental  Investigation  of  the  Base 
Pressure  Characteristics  of  Non  Lifting  Bodies  of 
Revolution  at  Mach  Numbers  from  2.73  to  4.98,  by  J.O. 
Reller  and  F.M.  Hanaker,  March  1955. 

20.  Hoerner,  S.F.,  Fluid-Dynamic  Drag,  Published  by  the 
Author,  1965. 

21.  Czysz,  P.,  Energy  Analysis  of  Propulsion  Systems  for 
High  Speed  Vehicles,  27th  Aerospace  Sciences  Meeting, 
January  9-12,  1989. 


49 


INITIAL  DISTRIBUTION  LIST 


No.  of  Copies 

1.  Defense  Technical  Information  Center  2 

Cameron  Station 

Alexandria,  Virginia  22304-6145 

2 .  Library,  Code  52  2 

Naval  Postgraduate  School 

Monterey,  California  93943-5002 

3.  Professor  Conrad  F.  Newberry,  Code  AA/NE  6 

Department  of  Aeronautics  and  Astronautics 

Naval  Postgraduate  School 
Monterey,  California  93943-5100 

4.  Commanding  Officer  2 

Attn:  LCDR  Benjamin  F.  Roper,  USN 

USS  John  Rodgers  (DD-983) 

FPO  AA  34092-1221 


50 


